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Abstract 

We  consider  the  solution  of  the  linear  systems  of  aigebraic  equations 
which  arise  from  elliptic  finite  element  problems  defined  on  composite 
meshes.  Such  problems  can  systematically  be  constructed  by  introducing 
a  bcLsic  finite  element  approximation  on  the  entire  region  and  then  belect- 
ing  subregions,  and  subregions  of  subregions  etc.,  where  the  finite  element 
model  is  further  refined  in  order  to  gain  higher  accuracy.  We  consider  con- 
jugate gradient  algorithms  where,  in  each  iteration,  problems  representing 
finite  element  models  on  the  original  region  and  the  subregions  prior  to 
further  refinement  are  solved.  We  can  therefore  use  solvers  for  problems 
with  uniform  or  relatively  uniform  mesh  sizes,  while  the  composite  mesh 
can  be  strongly  graded. 

The  central  theoretical  issue  concerning  these  so-called  iterative  re- 
finement algorithms  is  the  design  and  analysis  of  algorithms  for  which 
the  rate  of  convergence  is  independent  of  the  number  of  subproblems,  i.e. 
the  number  of  refinement  levels,  as  well  as  the  mesh  sizes.  Algorithms 
for  which  such  bounds  can  be  obtained  are  called  optimal.  In  this  paper, 
we  extend  our  previous  results  and  show  that  an  additional  important 
algorithm  is  optimal. 

As  in  our  previous  studies,  we  use  a  basic  mathematical  framework 
introduced  in  a  recent  study  of  a  variant  of  Schwarz'  alternating  algorithm. 
Some  tools  from  older  work  on  iterative  substructuring  methods  are  also 
used. 


1.  Introduction.  In  this  paper,  we  consider  the  solution  of  the  leirge  linear 
systems  of  algebraic  equations  which  arise  when  working  with  eUiptic  finite  element 
approximations  on  composite  meshes. 

Finite  element  models  on  composite  meshes  can  be  constructed  systematically 
in  a  framework  of  conforming  finite  elements;  cf.  Ciarlet  [3].  We  do  so  in  order 
to  be  able  to  use  a  number  of  technical  tools  which  aire  available  primarily  in  the 
conforming  case.  We  begin  by  introducing  a  basic  finite  element  approximation 
on  the  entire  region  and  we  then  repeatedly  select  subregions,  and  subregions  of 
subregions,  where  the  finite  element  model  is  further  refined.  The  resulting  Uneax 
system  is  solved  by  the  conjugate  gradient  method,  accelerating  the  convergence  by 
solving  so-called  standard  problems.  These  correspond  to  the  finite  element  models 
on  the  original  region,  prior  to  ciny  refinement,  and  those  on  the  subregions  prior  to 
further  refinement.  We  caji  thus  use  solvers  for  problems  with  uniform  or  relatively 
uniform  mesh  sizes,  while  the  composite  mesh  can  be  strongly  graded. 

This  approach  offers  a  number  of  advantages.  An  existing  code  can  be  up- 
graded, in  order  to  increase  the  accuracy  locally,  without  a  radical  redesign  of  the 
data  structures  etc.,  since  we  can  use  the  old  code  to  solve  one  or  several  of  the  stan- 
dard problems.  Issues  of  data  structures,  geometry  and  vectorization  are  generally 
easier  if  we  design  programs  for  composite  mesh  problems  in  terms  of  such  standard 
problems.  We  aJso  note  that  if  each  of  the  standard  problems  has  appreciatively 
fewer  degrees  of  freedom  than  the  composite  model,  then  we  might  benefit  from 
solving  a  number  of  smaller  problems  rather  than  one  large  one.  In  other  words, 
this  approach  can  also  be  viewed  as  a  divide-and-conquer  strategy. 

We  study  two  so-called  additive  algorithms.  These  methods  aie  particularly 
well  suited  for  parallel  computing  since  in  each  iteration  step  all  the  standard  prob- 
lems, on  the  the  different  levels  of  refinement,  can  be  solved  simultaneously.  One  or 
a  group  of  processors  czin  therefore  be  assigned,  in  a  straightforward  way,  to  each 
of  the  standard  problems.  Synchronization  between  the  processors  is  required  once 
in  each  iteration,  namely  when  the  residual  is  computed  and  assembled. 

The  systematic  study  of  iterative  refinement  methods  goes  back  at  least  to 
1983,  when  work  on  the  Fast  Adaptive  Composite  (FAC)  method  weis  begun  by 
Jim  Thomas  and  Steve  McCormick;  cf.  [14].  Issues  related  to  the  implementation 
of  this  method  on  paxedlel  computers  led  to  the  introduction  of  Asynchronous  FAC 
(AFAC)  methods  a  few  years  later;  cf.  Hart  and  McCormick  [10].  In  our  terminol- 
ogy, these  are  multiplicative  and  additive  algorithms,  respectively.  Numerical  ex- 
periments have  been  reported  for  model  cases  and  Richard  Ewing,  Steve  McCormick 
and  others  have  begun  to  test  the  algorithms  for  more  difiicult  problems  arising  in 
industry.  The  convergence  of  the  FAC  method  is  established  in  McCormick  et  al. 
[14],  [15]  under  a  certaiin  additional  regularity  assumption.  An  important  contribu- 
tion to  the  theory  is  given  by  Breimble,  Ewing,  Pasciak  and  Schatz  [2],  for  a  two 
level  FAC  algorithm.  Their  proof  of  optimality  is  regularity  free  and  so  are  ours. 
In  recent  papers  Bj0rstad  [1]  and  Mandel  and  McCormick  [12]  developed  a  theory 
for  two-level  multiplicative  and  additive  algorithms.  In  particular,  they  obtained 
interesting  results  concerning  the  relationship  between  the  spectra  of  the  two  meth- 
ods. In  another  recent  paper,  Mandel  and  McCormick  [13]  established  optimsJity 
for  a  multi-level  AFAC  algorithm  for  a  special  model  problem.  Recently,  Ewing, 
Lazarov  and  Vassilevski  [8]  obtained  interesting  bounds  for  finite  volume  schemes 


for  the  multi-level  multiplicative  algorithm.  In  our  own  earher  work  [19],  we  have 
established  optimality  for  this  zilgorithm  in  the  same  finite  element  framework  as 
used  in  this  paper.  The  main  result  of  this  paper  is  the  generalization  of  Mandel 
and  McCormick's  result  on  the  multi-level  case  to  general  conforming  finite  element 
models. 

Our  study  of  iterative  refinement  methods  begun,  at  the  Third  Copper  Moun- 
tain Conference,  with  the  discovery  that  the  FAC  and  AFAC  methods  have  a  struc- 
ture quite  similar  to  that  of  the  classical  Schwaxz  procedure,  see  Schwarz  [16],  and 
an  additive  variant  thereof  considered  in  Dryja  [4],  Dryja  and  Widlund  [5,6]  and 
Widlimd  [18].  That  work,  in  turn,  had  been  inspired  by  a  recent  paper  by  P.-L. 
Lions  [11]  in  which  a  variational  fraimework  for  the  classical,  multiplicative  Schwarz' 
method  is  developed  for  continuous  elliptic  problems.  Recently  we  have  also  shown 
that  some  so-called  iterative  substructuring  methods  also  can  be  analyzed  in  the 
same  basic  framework;  see  Dryja  and  Widlund  [6]. 

In  section  2,  we  introduce  the  finite  element  problems  on  composite  meshes, 
certain  projections  and  the  additive  algorithm  in  its  basic  form. 

In  Section  3,  we  consider  two  multilevel  additive  algorithms  and  provide  several 
bounds.  The  principal  result  is  that  one  of  these  methods  has  a  rate  of  convergence 
which  is  independent  of  the  number  of  refinement  levels  as  well  as  the  mesh  sizes. 

2.  Composite  Finite  Element  Problems  and  the  Basic  Iterative  Method. 

We  consider  lineax,  self  adjoint,  elliptic  problems  discretized  by  finite  element  meth- 
ods on  a  bounded  Lipschitz  region  Q  in  /?".  To  simplify  the  presentation,  we  assume 
that  the  differential  operator  is  the  Laplacian,  that  we  have  zero  Dirichlet  boundary 
conditions  and  that  we  use  continuous,  piecewise  linear  finite  elements.  Almost  all 
of  om-  results  can  be  extended  immediately  to  general  conforming  finite  element 
approximations  of  any  self  adjoint  elliptic  problem,  which  czm  be  formulated  as  a 
minimization  problem,  and  to  more  general  boundary  conditions.  The  continuous 
and  discrete  problems  axe  of  the  form 

a{u,v)  =  /(u),  y  V  eV, 

and 

aiuh,vh)  =  f{vH),  yvH  eV\  (1) 

respectively.  The  space  V  =  ifj(fi)  and  V''  C  V.  The  bilinear  form  is  defined  by 

a{u,  v)  =   I   Vu-Vv  dx  . 
Ja 

This  form  defines  a  semi-norm  |u|hi  =  {a{u,u)y/'^  in  H^{^)  and  a  norm  in  V. 
The  space  V^  is  defined  on  a  composite  triangulation,  which  might  be  the  result 
of  a  large  number  of  successive  refinements.  The  triangulation  of  Q  is  given  in  the 
following  way: 

We  first  introduce  a  relatively  coarse  triangulation  of  fi,  also  called  Qi,  and 
denote  the  corresponding  space  of  finite  element  functions  by  V''^ .  We  can  think 
of  this  space  as  having  a  relatively  uniform  (or  uniform)  mesh  size  hi.  Let  ^^2  be 
a  subregion  where  we  wish  to  increase  the  resolution.  We  do  so  by  subdividing  the 


elements  and  introducing  an  additionzd  finite  element  space  V''^.  We  assure  that 
the  resulting  composite  space  V'^^  +  V'''  is  conforming  by  having  the  functions 
of  V''*  vanish  on  5^2-  We  repeat  this  process  by  selecting  a  subregion  Q3  of  Q2 
cind  introducing  a  further  refinement  of  the  mesh  and  finite  element  space  etc..  We 
denote  the  resulting  nested  subregions  and  subspaces  by  17,  and  V'^'  respectively. 
Throughout,  we  assume  that 

and 

F''-'  nH^{n,)cV'''  CH^in,),      i  =  2,...,k. 

The  composite  finite  element  space,  on  the  repeatedly  refined  mesh,  is 

We  ctssume  that  all  the  elements  aire  shape  regulax  in  the  sense  that  there  is  a 
uniform  bound  on  h^/PK-  Here  /i/c  and  /?/<•  are  the  diameter  and  the  radius  of 
the  largest  inscribed  sphere  of  einy  element  K,  respectively.  Our  bounds  in  the 
theory  developed  in  this  paper  aiso  depend  on  the  shape  of  the  subregions  Hi.  Thus 
in  order  to  make  our  proofs  work,  we  cannot  allow  the  sets  Qi_i  \  Hi  to  become 
arbitrarily  thin  in  comparison  with  the  diameter  of  Hi-i.  We  edso  assume  that  the 
area  of  any  triangle  on  level  i  can  be  bounded  by  const.  q'~^ ,  times  the  area  of  the 
triangle  on  level  j  of  which  it  is  a  part.  Here  5  is  a  constant  <  1  and  j  <  i. 

The  finite  element  problem  is  defined  by  equation  (1)  and  the  corresponding 
stiffness  matrix  can  conveniently  be  computed  by  using  a  process  of  subassembly. 
Introducing  subscripts  to  indicate  the  domain  of  integration,  we  write 

a{u,v)  =  ani\nj(u,u)  +  an,\n3(u,u)  +  ...  +an»(u,u)  . 

The  stiffoess  matrices  corresponding  to  the  regions  Cl,  \  Q^+i ,  i  <  k  —  1,  and  Qk 
axe  computed  by  working  with  nodal  basis  functions  related  to  elements  relevTint 
to  the  subregion  in  question.  The  quadratic  form  corresponding  to  the  composite 
stiffness  matrix  is  then  the  sum  of  the  quadratic  forms  corresponding  to  Q,  \  fii+i, 
i  <  k  —  1,  and  Qk-  In  our  algorithms,  we  use  solvers  for  the  same  elUptic  problem 
on  the  subregions  Qi  ,  i  =  1, ...,  k,  and  the  meshes  corresponding  to  hj,  j  =  i  —  l,i  . 
The  corresponding  stiffness  matrices  cam  in  fact  be  obtained  at  a  modest  extra  cost 
during  the  assembly  of  the  composite  mesh  model.  When  we  refine  a  finite  element 
model  locally,  the  modified  stiffness  matrix  can  thus  be  obtained  by  replacing  the 
quadratic  form,  associated  with  the  subregion  in  question,  by  the  one  corresponding 
to  the  refined  model  on  the  szime  subregion.  It  is  therefore  relatively  easy  to  design 
a  method  which  systematically  generates  the  stiffness  matrices  for  all  the  standard 
problems  necessary  while  the  stiffness  matrix  of  the  composite  model  is  computed. 

The  fundamental  building  blocks  of  our  algorithms  are  the  P^,  j  =  i  —  l,i, 
the  projections  onto  the  spaces  V'^'  P|i7o(Q,).  We  note  that  if  j  =  t  —  1,  we  solve 
a  problem  on  Hi  with  a  coaxser  mesh  than  if  V'^'  were  used.  The  projection  P-  is 
defined  in  terms  of  the  unique  element  of  V'^'  fl  HQ{Cli),  which  satisfies 

a{Plvh,<l>h)  =  aivh,4>k),  V^^eF"'  nifoi(a).  (2) 


The  Schwarz  methods  can  be  defined,  in  general,  in  terms  of  a  set  of  subspaces 
V,  and  the  related  projections.  We  work  with  the  simplest  possible  polynomial  in 
the  projections  solving  the  equation 

PuH=(Pl+P2  +  ---  +  Pk)uh=9'H  (3) 

by  an  iterative  method.  Since  we  can  show  that  the  operator  P  is  symmetric,  with 
respect  to  the  bilinear  form,  and  positive  definite,  the  iterative  method  of  choice 
is  the  conjugate  gradient  method.  We  must  also  maJce  sure  that  equation  (3)  and 
equation  (1)  have  the  same  solution  by  choosing  the  correct  right  hand  side  for 
equation  (3).  Since  by  equation  (1) 

a(uh,(i>h)  =  fi<i>h)  , 

we  C2in  construct  the  right-hand  side  gj,  by  solving  this  equation  restricted  to  the 
different  subspaces  auid  by  adding  the  results.  It  is  similarly  possible  to  apply  the 
operator  P  of  equation  (3)  to  any  given  element  of  V''  by  applying  each  projection  P, 
once.  Most  of  the  work  ,  in  particular  that  which  involves  the  individual  projections, 
can  be  ccirried  out  in  parallel. 

It  is  well  known  that  the  number  of  steps  required  to  decrease  an  appropriate 
norm  of  the  error  of  a  conjugate  gradient  iteration  by  a  fixed  factor  is  proportioned 
to  y/n,  where  k  is  the  condition  number  of  P;  see  e.g.  Golub  and  Van  Loan  [9].  We 
therefore  need  to  establish  that  the  operator  P  of  equation  (3)  is  not  only  invert ible 
but  that  satisfactory  upper  and  lower  bounds  on  its  eigenvalues  cam  be  obtained. 

3.    Two  Additive  Algorithms.   An  additive  algorithm  is  defined  by  specifying 
the  subspaces  Vj  or  alternatively  the  projections  Pi.     Before  we  discuss  specific 
algorithms,  we  maice  some  remarks  on  estimating  the  eigenvalues  of  P  from  above. 
The  upper  bound  on  the  spectrum  is  obtained  by  bounding 

aiPvh,Vh)  =  a{PiVh,VH)  +  a{P2Vh,Vh)  +  ---  +  a{PkVh,Vh)  ,      . 

from  above  in  terms  of  a{vh,Vh).  We  can  use  Schwarz'  inequality  and  the  fact  that 
P,  is  a  projection  to  prove  that  each  term  is  bounded  by  a{vh,Vh)  and  thus  the 
spectrum  of  P  is  bounded  from  above  by  k.  Our  goal,  however,  is  to  establish 
a  uniform  bound  on  the  condition  number.  This  can  be  done  if  the  terms  are 
orthogonal  or  almost  orthogonal;  cf.  Theorem  1. 

We  consider  two  different  algorithms  and  distinguish  between  them  by  using 
the  superscript  1  and  2.  The  most  natural  algorithm  uses  the  projections  P^^^  =  Pi, 
where  the  projection  operators  P/  are  defined  in  equation  (2).  The  condition  number 
of  this  algorithm  grows  linearly  with  k.  This  can  be  shown  as  follows.  We  have 
already  remarked  that  the  eigenvjilues  of  P  always  are  bounded  from  above  by  k. 
This  bound  is  attained  if  V'**  D  H^(^k)  is  not  empty,  i.e.  when  the  mesh  size  is  fine 
enough.  Any  such  function  belongs  toV''',i  =  l,2,...,k,  and  is  exactly  reproduced 
by  each  of  the  projection  operators.  It  is  thus  an  eigenfunction  with  eigenvalue  k. 
Similarly,  any  function  which  belongs  to  F'"  nH^{^i\Q,2)  is  an  eigenfunction  with 
eigenvalue  1.    By  using  the  same  construction  as  in  the  proof  of  Theorem  2,  we 


can  show  that  the  eigenvalues  of  P^^^  are  bounded  from  below  by  a  constant.  The 
condition  number  of  P^^'  is  therefore  of  order  k. 

We  are  principally  interested  in  the  method  defined  by  P^  '   =  P,'  —  P^^j   , 

i  <  k  —  I,  and  Pf.  =  Pj^  as  the  basic  projections  in  equation  (3).  It  easy  to  show 
that  these  operators  axe  projections  and  we  will  see,  in  the  proof  of  Theorem  2, 
that  the  composite  finite  element  space  V^  is  the  direct  sum  of  the  corresponding 
subspaces  ¥■ 

Our  main  new  result  is  the  lower  bound  given  in  Theorem  2.  For  completeness, 
we  also  include  our  older  proof  of  the  upper  bound  on  the  spectrum  of  P. 

Theorem  1.  The  eigenvalues  of  P^^^  are  uniformly  bounded  from  above  by  a 
constant. 

Before  we  turn  to  the  proof  proper,  we  observe  that  the  unbalance  of  the  first 
method,  which  resulted  in  the  ajnplification  of  certaiin  functions  by  a  factor  A:,  is  no 
longer  possible.  By  regrouping  terms,  we  also  see  that 

P(2)  =  p^  +  (p2  _  p^l  )  +  .  .  .  +  (^pk  _  pfc-1  )   _  (4)  . 

All  the  terms,  except  the  first,  are  similar  to  multigrid  corrections  since  they  each 
represent  the  difference  between  two  solutions  on  the  same  subregion,  using  two 
different  mesh  sizes. 

By  using  the  representation  of  P^^\  given  in  equation  (4),  we  obtain 

P^^'>Uh  =  PIUH  +  (P|  -PI)UH  +  ---  +  {Pt  -  P^-')UK 
=  Ui  +U2  +  •■■  +  Uk- 

We  show  that  these  terms  axe  increeisingly  orthogonal: 

a{ut,Um)  <  const.  gi'"'"'(a(u/,Uf))^/^(a(um,Um))^''^, 

where  ^i  <  1,  uniformly.  This  is  a  so-called  strengthened  Cauchy  inequality.  With- 
out loss  of  generality,  we  assume  that  i  <  m.  It  is  easy  to  show  that  a{um,Vh)  =  0, 
Vu/i  G  y''"'-ini?o(nm),  i-e.  u„j  is /im-i— harmonic  on  fim-  Let  u^  =  ^  a,(^,'"  where 
the  (f>- '  are  nodal  basis  functions  in  V''' .  Any  such  basis  function  is  orthogonal  to 
Um,  if  its  support  does  not  intersect  dQrni  since  then  either  its  support  belongs  to 
^e  \  ^m  or  (jf»|"  €  F''""-!  n  H^{Qrn)-  Thus  the  only  contributions  to  a(u„i,  u^)  come 
from  the  triangles  of  V'''  which  intersect  dQm-  Consider  one  such  triangle  T.  The 
restriction  of  ug  to  T  is  a  linear  function  and  thus  has  a  constant  gradient.  It  can 
also  be  written  as  a  hnear  combination  of  basis  functions  in  F'''"-'.  The  support 
of  most  of  these  do  not  intersect  dQ.m  and  therefore 

aT{u,n,ut)  =  aTnsium,ut)  <  (aT("m,Um))^^^(arns(«^u^))^''^ 
<  const.  9l'""''ar(um,m)^/^ar(u/,u^)^/^  . 

Here  S  denotes  the  union  of  the  small  triangles  that  are  next  to  dQm-  We  also  use 
our  assumptions  on  the  triangulations,  given  in  section  2,  and  the  fact  that  the  gra- 
dient of  ut  is  constant  on  T  and  that  therefore  its  contribution  to  the  quadratic  form 


is  directly  proportional  to  the  axea  of  integration.   The  proof  of  the  strengthened 
Cauchy  inequality  is  completed  by  summing  over  the  triangles  of  V^' . 
We  now  note  that 

it 


Here 


A  =  const. 


/I  91  ll 

qi        1        Qi 

\  -,*-!  .,*-2  ^t-3 

^9i  Q\  1i 


and 


It  is  easy  to  show  that  the  EucHdecin  norm  of  .4  is  uniformly  bounded.  Thus 

k  k 

a(P<2'u/,, P^^^Uft)  <  \A\e^  ^  a{ui, u.)  <  const.   ^  a{u„u,). 

;=1  «=1 

To  complete  the  proof,  we  note  that 

a(ui,u,)  =  aiPiuh  -  Pj-'^uk,  P!uh  -  P/"^uO 
=  a(P/u,,-P;-iuA,P/uA), 

since  P^Uh  —  P*~^Uh  is  /ii_i— harmonic.     By  using  elementary  properties  of  the 
projections,  we  find  that  a(ui,Ui)  =  a(uj,u/,).  Therefore 

it  k 

and  thus, 

a(p(2)y^  p(2)^^^  ^  ^^^^^   a(p(2)«h,Uh), 

from  which  the  upper  bound  on  P^^^  follows. 

We  now  complete  our  proof  of  the  optimality  of  the  second  additive  algorithm 
by  proving  a  uniform  lower  bound. 

Theorem  2.   The  eigenvalues  of  P^"^^  are  uniformly  bounded  from  below. 

We  begin  the  proof  by  establishing  some  auxiliary  results. 
A  lemma,  essentially  given  in  Lions  [11],  has  proven  very  useful  for  obtaining 
lower  bounds.  Since  the  proof  is  quite  short,  it  is  included. 

Lemma  1.  Letuh  =  J2i=-[  Uh,i,' where  u/,,,  GV,,  be  a  partition  of  an  element  of 
V"  =  Vi+...  +  Vjt.  If  the  representation  can  be  chosen  such  </iai  Xl.*=i  (^iuh,i,uh.i)  < 
Co2a(u/„Uft),  VuA  G  F^  t/ien  A^.„(P)  >  Cq-^. 
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Proof:  By  elementary  properties  of  symmetric  projections  and  the  represen- 
tation of  u/i  as  a  sum,  we  find  that 

k  k  k 

a{uh,Uh)  =  ^a{uh,Uh,i)  =  '^a{uh,P,Uh,i)  -  '^a{PiUh,Uh,i)  . 

.=  1  •=!  i=l 

Therefore, 

k  k 

aiuh,UH)  <  {J2aiP,UH,P,UH)Y^H^a{uh,„UH,,)y/'^  . 
1=1  .=1 

By  the  assumption  of  the  lemma 

k  k 

i=\  1=1 

and  the  lemma  is  established. 

In  addition  to  the  projections  P/,  we  also  use  another  fcimily  of  operators, 
ifj^j,  defined  by 

a{Hi^,Vh,WH)  =  0,    Vu;,.  e  V'"^  n  H^i^j+i). 

It  is  natural  to  czdl  Hj^^Vh  the  /ij— harmonic  extension  of  Vh.  to  f2j+i,  since  it 
is  the  solution  in  V'^'  of  a  discrete  Dirichlet  problem  with  zero  right  hcmd  side  and 
with  boundary  data  on  dQj+i  given  by  Vh- 

The  following  lemma  is  a  consequence  of  the  extension  theorem  given  in  Wid- 
lund  [17]. 

Lemma  2.  There  exists  a  constant  C  ,  which  is  independent  of  Uh  and  the 
mesh  sizes,  such  that 

an,+,{Hj^iUh, Hj^iUh)  <  CaQ.\Q._^_^{uh,Uh),  Vufc  €  V''. 

We  will  not  discuss  the  proof  of  this  lemma.  We  note  that  the  constant  C 
necessarily  blows  up  if  we  let  the  area  of  Q.j  \  Q.j+\  shrink  to  zero,  keeping  Qj+ij 
fixed.  Such  situations  are  not  of  any  real  interest  in  our  application. 

Lemma  3.  If  the  restriction  of  Vh  to  fij+i  belongs  to  V^' ,  then 

Vh  =  Hj^^Vh  +Pj+iVk 
and  the  two  terms  are  orthogonal  in  the  sense  of  the  bilinear  form. 
Proof:  It  is  easy  to  see  that 

a((/  -  ifj+,K,<A/.)  =  a{Pj^,VH,<l>h),  V<^fc  €  V"^  D  H^iQj+i)  . 
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Since  it  follows  from  the  assumption  that  {I  —  Hj^i)vh  €  V*'  n^^(Qj+i ),  we  obtain 

Proof  of  Theorem  2:  To  prepare  for  the  use  of  Lemma  1,  we  partition  an 
arbitrary  u/,  €  V^ 

k 

i=l 

where  V^^     is  the  range  of  P^    .  The  terms  u/,,^  are  given  by 

ut,^,  =  Hl^i(uh  -  S,-iUh)  ,      i<k-l,  (5a) 

where 

j 

SjUh  =  2_^  "ft.' ' 
1=1 

and  by 

Uh,k  =  Uh  -  Sk-iUh.  (56) 

The  easiest  way  to  see  that  the  sum  of  the  Uh,i  indeed  equals  u/,  is  to  consider 
the  values  on  the  sets  Q,,  \  Q,i+i  separately.    To  prove  that  if'^j(u/,  —  Si-iUh)  G 

.  y,  ,  we  apply  the  operator  P/  —  P/^j  to  this  function.  Since  this  function  is 
/li-harmonic  on  Qi+i  ,  P,*.^.iU/,,,  =  0  .  The  argument  is  completed  by  observing  that 
■H''^.i(u/,  —  Si-iUh)  G  V**  .  Therefore  (P/  —  P'^.i)u/,,i  =  u/,,,  and  the  correctness  of 
formula  (5)  is  established. 

The  partial  sums  SjU^  can  also  be  written  as 

SjUh  =  Sj-iUh  +  Hj^iiuh  -  Sj-iUh)  =  Hj^iUh  +  (/  -  H^j^i)Sj-iUh. 
Therefore  by  Lemma  3, 

SjUh=Hj^,Uh  +  P^+,Sj.iUH,  (6) 

where  the  two  terms  axe  orthogonal. 

We  now  proceed  to  estimate  the  norm  of  the  different  terms  in  the  represen- 
tation of  Uh  given  by  equation  (5).  We  first  consider  Uh,k  G  H^{Q.k),  using  the 
orthogonality  of  the  two  terms  of  equation  (6), 

aiuh,k,uh,k)  =aQ^{uh,k,Uh,k) 

<  2an»(uft,ufc)  +  2aQ^{Sk-iUh,Sk-iUh) 

=  2an,{uh,Uh)  +  2an,{H^-'uh,H^-'uh) 
+  2an,(Pfc*-'Sfc_2tx/.,P*-^5fc_2u„)  . 

The  next  term  can  be  estimated  as  follows: 

a(uh,k-i,Uh,k-\)  =  au^_^{uh,k-uUh,k-i) 

<^<i^..AH'k-'uh,H',-'uh) 

+  ^(lQ.-dH'k-'Sk-2Uh,H',-'Sk-2Uh). 
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By  the  definition  of  H^~\ 

Thus 
a{uh,k,Uh,k)  +  a{uh,k-i,Uh,k-i) 

+  2an,iPt'Sk-2UH,Pl:-'Sk-2Uh)  +  2an,.AHk~'Sk-2nH,H'k-'Sk-2Uh). 

By  Lemma  3,  the  sum  of  the  last  two  terms  is  equed  to  2an^_^{Sk-2Uh,Sk-2Uh)- 
The  same  arguments  can  now  be  repeated  and  we  obtain 

a{uh,k,Uh,k)  +  a{uh,k-i,Uh,k-i)  +  a{uk,k-2,Uh,k-2) 

<2an,_,{uh,Uh)+AaQ,{H^^-'^Uh,H^-'^Uh) 
+  4anfc_i(-H'iZ^w/,,  H^Zi^h)  +  2aQ^_j{Sk-3Uh,  Sk-jUh) 

and  finally 

a{uh,k,Uh,k)  -\ \-ci{uh,i^^h,i) 

By  using  Lemma  2  and  by  combining  the  contributions  from  the  regions  Q;-:  \ 
Qj,j  =  2,--,k,  into  one,  we  see  that  Cq  <  2  +  4C  .  Theorem  2  now  follows  directly 
from  Lemma  1  . 
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